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Abstract 

The efficient realization of linear space-variant (non-convolution) 
filters is a challenging computational problem in image processing. In 
this paper, we demonstrate that it is possible to filter an image with a 
Gaussian-like elliptic window of varying size, elongation and orienta- 
tion using a fixed number of computations per pixel. The associated 
algorithm, which is based on a family of smooth compactly supported 
piecewise polynomials, the radially-uniform box splines, is reaUzed using 
pre-integration and local finite-differences. 

The radially-uniform box splines are constructed through the re- 
peated convolution of a fixed number of box distributions, which have 
been suitably scaled and distributed radially in an uniform fashion. The 
attractive features of these box splines are their asymptotic behavior, 
their simple covariance structure, and their quasi-separability. They 
converge to Gaussians with the increase of their order, and are used 
to approximate anisotropic Gaussians of varying covariance simply by 
controlling the scales of the constituent box distributions. Based on the 
second feature, we develop a technique for continuously controlling 
the size, elongation and orientation of these Gaussian-like functions. 
Finally, the quasi-separable structure, along with a certain scaling prop- 
erty of box distributions, is used to efficiently realize the associated 
space-variant elliptical filtering, which requires 0(1) computations per 
pixel irrespective of the shape and size of the filter 

Index - Space-variant filter, Finite-differences, Running-sums, Anisotropic 
Gaussian, Box spline, Zwart-Powell (ZP) element. 

1 INTRODUCTION 

The most widely used smoothing operator in image processing is the Gaus- 
sian filter. As far as isotropic Gaussians are concerned, a fast implementation 
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is achievable simply by decomposing the filter into two orthogonal 1-D Gaus- 
sians operating along the image axes. The 1-D filters are in turn implemented 
using efficient recursive algorithms, e.g., the ones proposed by Deriche [5] 
and Young et al. [20]. We refer the readers to this survey article [16] for an 
exhaustive account of such recursive schemes. 

A fundamental limitation of isotropic filtering is that it does not take into 
account the anisotropic nature of image features, which results in blurring of 
oriented patterns and textures. The development of fast anisotropic filtering 
in general, and anisotropic Gaussian filtering in particular, have therefore 
gained momentum over the past decade. Worth mentioning in this regard 
is the work of Geusebroek et al. [7], who developed an efficient recursive 
technique based on the factorization of the 2D -anisotropic Gaussian into two 
1-D Gaussians, one along the image axes and the other along a generally off- 
grid direction. A drawback of this technique is that one has to interpolate the 
image along the off-grid direction to implement the corresponding 1-D filter. 
To avoid interpolation and, in effect, to improve the spatial homogeneity and 
the Gaussian-like structure of the filters in [7], Lam et al. came up with the 
alternative "triple-axis" solution. Instead of using two directions, they chose 
to decompose the anisotropic Gaussian into three 1-D Gaussians operating 
along one ofthe four cardinal directions: the horizontal, the vertical, and 
the two diagonals [10]. The focus of these papers has largely been on 
space-invariant filtering, where the entire image is convolved with a single 
anisotropic Gaussian. On the other hand, a variety of space-variant filtering 
strategies have also been developed, including image statistics driven filtering 
[11], non-linear diffusion filtering [14, 19] and gradient inverse-weighted 
filtering [18], to name a few. 

1.1 Linear space- variant filtering 

In this paper, we focus on the paradigm of linear space-variant filtering using 
Gaussian-like kernels of different shapes and sizes. From a purely discrete^ 
perspective, this calls for the design of a family of Gaussian filters {^Ai^llA^ 
so that, given an input image f[n], one can evaluate the filtered samples 
f[n\ through the averaging 

/>]=E/[%A(n)[n-fc]. (1) 
k 

^We associate the term "discrete" with functions defined on the Cartesian lattice Z'^, where 
d is the dimensionahty of the function. 
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The parameter A(n), which specifies the covariance of the filter appHed 
at location n, allows one to continuously adjust the scale, orientation and 
elongation of the filter in keeping with the anisotropy of the local image fea- 
tures. There are however certain practical problems involved in an efficient 
realization of (1). It is obvious that (1) cannot be written as a convolution, 
and hence cannot be realized using the FFT algorithm. In fact, the available 
options are either (i) to compute the filters g\ [n] by sampling the contin- 
uous Gaussian on-the-fly, or (ii) to discretize A a priori, and to store the 
pre-compute filters in a look-up table. The problem with the former is that 
it proves to be extremely slow for wide kernels, while the latter restricts 
the control on the anisotropy of the filters. By appropriately modifying the 
recursive filtering strategies in [5, 20], Tan et al. developed an algorithm for 
reahzing (1) for the particular case where {^a W}a are isotropic [16]. 

Spline kernels can also jrield efficient algorithms for space-variant filter- 
ing, particularly when the space-variance is in terms of the scale (or size) 
of the kernel. For instance, Heckbert proposed an algorithm for adaptive 
image blurring using tensor-products of polynomial splines, where the image 
is filtered with kernels of various scales using repeated integrations and 
finite-differences [8] . Based on similar principles, namely, the scaling prop- 
erties of B -splines, Muiioz-Barrutia et al. have developed an algorithm for 
fast computation of the continuous wavelet transform of 1-D signals [13]. 
Recently, the method was extended to perform space-variant filtering using 
Gaussian-like functions of arbitrary size, which can be elongated along the 
image axes [12]. To achieve this, the authors choose to approximate the 
Gaussian using separable B -splines. We propose to take this approach one 
step further. In particular, we overcome the limited steerability and ellipticity 
of the separable B -splines by considering certain quasi-separable analogues 
of B-splines, the so-called box splines [2]. As was demonstrated for the sepa- 
rable filters in [12], we show that these quasi-separable box splines can also 
be used to approximate the Gaussian, and that the associated space-variant 
filter can be decomposed into recursive pre-filters and scale-dependent finite 
difference filters. These together allow us to achieve a fast space-variant 
filtering of images using elliptical Gaussian-like filters. 

To date, there have only been few applications of such multivariate splines 
in the fields of image processing and computer graphics. Noteworthy among 
them are the works of Richter [15] and Asahi et al. [1], concerning the 
development of image approximation and reconstruction algorithms, and 
that of Condat et al. [4] and Entezari et al. [6], concerning the development 
of interpolation formulas for hexagonal and BCC lattices. 
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1.2 Main idea 



In this contribution, we propose a fast space-variant filtering algorithm using 
a family of Gaussian-like box splines whose size, elongation and orientation 
can be continuously controlled. The attractive feature of our approach is 
that we use a continuous-discrete formalism which avoids the necessity of 
sampling a continuously-defined filter on-the-fly, or of storing a discrete set of 
pre-computed kernels. The developments in the paper are centered around 
two main ideas, which are as follows: 

(1) The use of quasi-separable box splines. The construction of bivariate 
box splines, conceived as the "shadow" of TV-dimensional (N > 2) polj^opes 
in 2-dimensions, often turns out to be rather intricate (see [2] for instance). 
In this paper, we consider an alternative straightforward recipe for con- 
structing box splines, namely, through repeated convolutions of dilated and 
rotated box distributions (see Fig. 3). In particular, we realize the so-called 
radially-uniform box spline P^{x) through the convolution of N rotated box 
distributions, where a = (04, . . . , a^v) is a scale-vector with au being the scale 
of the box distribution positioned along the direction {k — l)Tr/N (see §2.2 
for a precise definition). 

The reason why the radially-uniform box splines are of interest in the 
current context is twofold. The first of these is that we can make them 
arbitrarily close to a Gaussian by increasing (see §2.2.3). The second 
reason, which has a more practical significance, is that we can continuously 
control their size, elongation, and orientation simply by acting on the scales. 

(2) An efficient strategy for space-variant averaging. To convey this idea, 
we examine the following formula 



for computing the space-variant averages of a discrete 1-dimensional signal 
f[n\. We interpret the factor W{n), which controls the amount of smoothing, 
as the size of the "discrete box filter" applied at location n. The disadvantage 
of using (2) is that its computational cost scales linearly with W{n), which 
even gets worse in higher dimensions. This can be circumvented (with a 
mild interpolation cost) by considering instead the space-variant averaging 




W{n) 



(2) 




(3) 
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n-a{n)/2 n + a{n)/2 



Figure 1: Computation of the space-variant average f[n]: The signal /(x) is 
first localized (hatched zone) using the shifted box function fia{n) ~ ^rid 
then the area is computed. The central idea of our algorithm is to determine 
this area by taking the finite-difference of the primitive of f{x). 

where we have replaced the discrete signal f[n] by its interpolated version 
f(x), and the discrete box filter by the box function 




l/a for -a/2 < x < a/2, 
otherwise. 



The main advantage of this formulation is that we can realize (3) using 0(1) 
computations per position, independent of the size of a(n). This is based on 
the observation that (3) can computed by first evaluating the primitive 

F{x) = r f{y)dy, 

J — oo 

and then using the formula 

/N = f^(" + ain)/2) - F{n - a(n)/2)) , (4) 

which requires one addition and multiplication per position. This idea is 
illustrated in figure 1. The other advantage is that, as opposed to the integer- 
valued window W{n) in (2), this gives access to the real-valued scale a(n) 
for continuously controlling the smoothing. Indeed, if f{x) is integrable (at 
least locally), then it can be shown that the use of small scales results in less 
smoothing, namely that f[n] — y f[n] as a(n) — y 0, whereas f[n] can be 
made negligibly small by making a(n) sufficiently large. 
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The contribution of this paper is the generalization of the filtering strategy 
in (3) to the bivariate setting using the radially-uniform box splines. In 
particular, given a discrete image f[n\, we consider the space-variant filtering 



where f{x) represents a suitable interpolation of /[n]. The significance of the 
quasi-separable characterization of j3^{x) in terms of the box distributions 
is that it allows us to relate (5) to the 1-D problem in (3). Indeed, we 
demonstrate in §2.2 that (5) can be implemented using an appropriate 
bivariate extension of pre-integrations and finite-differences, together with 
few evaluations of a fixed piecewise polynomial function (the coefficients 
of which are pre-computed). Although the derivation of the algorithm is 
rather involved, the final solution turns out to be remarkably simple (see §3, 
Algorithm 1), and easy to implement. 

1.3 Notations 

We use f{u)) to denote the Fourier transform of a function f{x) on R*^, 
specified by f{uj) = f{x) exp {—juj^x)dx. We use /(• — s) to denote 
the function obtained by translating f{x) by s. The convolution of two 
functions /(a?) and ^(a;) is given by (/ * g){x) = J^a f{s)g{x — s)ds. The 
notation ®^^ifk{x) is used to denote the convolution of a collection of 
functions fi{x), . . . , /at (a;); the order of the convolutions is immaterial. We 
suppress the domain of an integral (or summation) if this is obvious from 
the context. For a bivariate function f{x) = f{xi,X2), the partial derivative 
along Xi is denoted by dif{x). Given operators Ti and T2 on a domain 
V, we use Ti o T2 (often simply T1T2] to denote their composition: (Ti o 
T2){f) = Ti(T2(/)) for every / in V. We use M" to denote the (n - 1)- 
fold matrix multiplication of M with itself. The integral J'M.{x)f{x)dx, 
corresponding to a real-valued function f{x) and a matrix- valued function 
M(a;) on R^, denotes a matrix of the same dimension as M(a;), whose 
(i, j)-th component is given by J M.i^j{x)f{x)dx. If P and Q are constant 
matrices, we then have / PM(a;)Q/(a;)da; = P(/ M(a;)/(a;)da;)Q. The 
notation f{x) = 0{g{x)),x — > 0, signifies that there exists a constant C 
(independent of x) such that |/(x)| < Cg{x) for sufficiently small x. The 
space of bivariate finite-energy signals is denoted by L-^(R^), or simply as L^; 
it is equipped with the norm ||/||l2= [/ |/(a;)p(ia;]^/^. The Dirac distribution 
is denoted by 5{x). 




(5) 
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Ar'/3i(^) 



= A„A-i/3i(x + r) 



Figure 2: Box function rescaling through "addition and subtraction" of the 
unit box function: The step function is first reproduced from the unit box 
using the running-sum, and then the appropriate finite-difference is appHed 
to recover the rescaled box function. 

2 SPACE-VARIANT AVERAGING 

We now derive (4) using an operator-based formalism. This helps set up the 
framework needed for the subsequent generalization of the idea to higher 
dimensions and multiple orders in §2.2.2. 

2.1 Realization of (3) 

We consider the finite-dijference (FD) and the running-sum (RS) of a function 
f{x), which are respectively defined by 

Aa/(x) = ^(/(x)-/(x-a)), (6) 

and 

oo 

A-'f{x) = bY,f{x-bk). (7) 

The positive real numbers a and b are the scales of the operators^. We note 
that the operators and Aj^^, which takes f{x) into Aaf{x) and A^^f{x) 
respectively, are linear and translation-invariant, and that when b is an 
integer, Aj^^ can also be applied to a sequence f[n] through the well-defined 
operation Aj^^/[n] = 6^^q f[n — bk]. In particular, g[n] = A^-^fln] can be 
implemented efficiently using the simple recursion g[n] = g[n — 1] + f[n], 
under appropriate boundary conditions [13]. 

^The notation AjJ"^ is justified by the fact that Aa A^^ acts as the identity operator when 
a = b. 
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The significance of these operators is that we can relate the variable-size 
box functions in (3) to the unit-width box function using the transformation 

/3a(x) = AaA^Vi(2; + r) (8) 

where r = (a — l)/2. In particular, this means that box functions of variable 
widths can be derived from a fixed box function through the successive 
applications of running-sums and finite-differences (see Fig. 2). To derive 
(8), we note that A^^/3i{x) = Xlfclo - k) = u{x + 1/2), where the step 
function u{x) equals 1 for x > and zero otherwise. The desired result 
follows immediately: 

AaA5-^/3i(x + r) = ^ {u{x + a/2) - u{x - a/2)) = Pa{x). 

We use (8) to derive the algorithm for computing (3) as follows: Fix 
an arbitrary position n and the corresponding a(n) in (3), and consider the 
function 

six) = J f{y)Pi{x - y)dy. 

We claim that /[n] = A(,(^)Aj^^s(n + r). Indeed, following the linearity and 
translation-invariance of Aa{n)^b^> using (8), we can write 

A„(„)A]-is(x + r) = j f{y)[l^a{n)^i^Pi{x + T-y)]dy 
= j f{y)Pain){x -y)dy, 

which establishes our claim. 

Now if the input signal is discrete, of the form f{x) = Xlnez /N"^!^ ~ 
n), then s{x) can simply be written as s{x) = ^ f[n](3i{x — n). A simple 
manipulation then shows that A^^s(a;) = Y^g[n\l3i{x — n), where g[n\ is the 
running-sum of f[n]. Thus, denoting the piecewise-constant interpolation of 
g[n] by F{x), we obtain 

/[n] =A„(„)F(n + r). 
This leads us to the following two-step algorithm for realizing (3): 

(1) (Space-invariant step) Compute g[n] = A'^^f[n] using the recursion 

g[n] = g[n - 1] + f[n]; 
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(2) (Space-variant step) For every position n, set f[n] = A„(„)F(n + r), 
where r = (a(n) — l)/2. 

The steps of the algorithm can be visuaHzed for the particular case when 
the input is an impulse and when a(n) = a for every n using Fig. 2. The 
second and third plots in this figure then correspond to steps (1) and (2) of 
the algorithm, respectively 

The remarkable fact about the algorithm is that the space-variant aspect 
of the transformation f[n\ i-> f[n\ gets transferred to the scale-dependent 
operator Aq, which in turn is implemented at a fixed computational cost per 
pixel, namely, one addition and multiplication per position. We would also 
like to note that (8) can more generally be expressed as 

/3aix) = AaA-'/5,{x + T) (t = (a - 6)/2). (9) 

The significance of this relation is that, if the lattice spacing b is different 
from unity, one can still realize the running-sum (without interpolation) by 
replacing the operator A^^ by A^^. We will use this in the sequel. 

2.2 Bivariate extension 

2.2.1 Radially-uniform box splines 

We now extend the space-variant filtering strategy discussed in the previous 
section to the bivariate setting, where the additional notion of directionality is 
appropriately addressed. As a first step, we devise an appropriate directional 
extension of the box function. In particular, corresponding to a real-valued 
scale a and direction 0, we define the directional box distribution ipa,e{x) as 
the tensor product of the box function f3a{x) and the Dirac distribution 5{x) 
operating along orthogonal directions. 

Here the orthogonal directions are specified by the unit vectors ue = 
(cos 0, sin 0), and ug^ = {— sin 6 , cos 6) . The scale a controls the amount 
of smoothing applied along the orientation of the box distribution, whereas 
no smoothing is applied along the transverse direction. The idea then is 
to construct the box spline by convolving an arbitrary number of such di- 
rectional box distributions (cf. Fig 3). Thus, corresponding to an integer 
iV > 1, a set of real-valued scales oi, . . . , a^v, and uniform rotation-angles 
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Figure 3: Construction of the radially-uniform box spline through the convo- 
lution of four directional box distributions. (A) The four box distributions, 
distributed uniformly over [0, vr), were assigned equal scales in this example; 
(B) Scan profile along 6 = ir/S. 

= {k — 1)71 /N, k = 1,. . . ,N, we specify the radially-uniform box spline 
through the convolution 

/3f (x) = [ifa^fi^ *■■■* ^Paj^,ej^)ix). (10) 

We shall refer to N and the tuple a = (ai, . . . , oat) as the directional-order 
and the scale-vector of the box spline, respectively. 

Following definition (10), it can be verified that (x) is a piecewise 
polynomial of degree < — 2, where the partitions are specified by lines 
running along the directions 6*1, ... , 9^. Moreover, I3^{x) is symmetric with 
respect to the origin, and is compactly supported on a convex A^-sided 
polygon consisting of the points 

N 

{^tkakue, : -1/2 <tk < I/2}. 

k=l 

The radially-uniform box splines are non-separable for N > 2. However, 
in keeping with the spirit of the underlying tensor construction, the term 
quasi-separable would be more appropriate. The scale-vector a plays a vital 
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role in determining the size and shape of the box spHne. It is clear that 
the box spline can be arbitrarily elongated along the principal directions 9k 
{I < k < N) simply by rescaling the box distribution (fa^Mk' that is by making 
Qk large compared to the other scales. Moreover, we will demonstrate in the 
sequel that one can elongate the box spline along any arbitrary direction by 
appropriately acting on the scale-vector The role of the directional-order is 
more subtle; it determines the degrees of freedom available for controlling 
the geometry of the box spline and also its regularity (smoothness) . We will 
discuss these aspects in detail for the particular four-directional box sphne 
(AT = 4) in §2.2.3. 



2.2.2 Realization of (5) 

We now formulate the algorithm for realizing (5) by appropriately extending 
the domain of definition of the FD and the RS operator to bivariate functions. 
The main idea is to derive a relation similar to (8) for the radially-uniform 
box splines. Thus, corresponding to positive real-valued scales a and b, and 
direction < ^ < tt, we consider the directional finite-difference 

Aa,ef{x) = ^(^f{x)-f{x-au0)), (11) 
and the directional running-sum 

oo 

^b,lfi^) = bY,fi^- kbue). (12) 

k=0 

In keeping with the definition of the box spline, the radially-uniform FD and 
RS operators, and A^^, are then specified by the combined action of 
(11) and (12) along the directions 6k = {k — 1)it/N. In particular, we set 

A^ = A„,,e, o...oA„^,e;v> (13) 

and 

Ab'^ = A,-;,^o...oA-j_,^, (14) 

where the scale-vectors a = (ai, . . . , cat) and b = {hi, ... ^h^) specify the 
scale along each direction. The operators A^ and A^^ are closely related to 
the radially-uniform box splines. Indeed, it can readily be verifed that 

A^ A-^ = K,,eAble, ° " " " ° ^a^^fir^^llfi,, 
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and that using (9) we can write 



^a,e\l<fb,e{x + TUe) = (Pa,e{x). 

Thus, if we let r = X] TkUe^, where = (a^ — 6fe)/2, then following defini- 
tions (10), (13), and (14), we see that 

,N \-N oN 



= ®k=i^a,,e,A^^]o^ipb^,edx + Tkue^) 



k=i 'Pbkfikix + Tkuej 



This provides the following crucial connection between the box splines and 
the 2-D operators. 

Proposition 2,1. The box spline P^ix) can be expressed as 

/3^{x) = A^A-''/3^{x + t). (15) 

Before discussing the filtering algorithm, we briefly elaborate on the 
implementation of and A,^^. We can show that (13) can be written as 

Kf{x)= E mf{x-xi), (16) 

i=0 

where Wi = (-l)'?i+-+9^(ai • • • oat)-^ and Xi = J2k=i IkakUe^^, and the 
index i taking values between and (2^ — 1) is the decimal counterpart of 
the binary number (qn, qn-i, ■ ■ ■ , qi), which takes values from (0, . . . , 0) to 
(1, . . . , 1). We identify the points Xi with the vertices of an affine mesh, and 
Wi with the corresponding mesh taps (cf. Fig. 5). 

As far as the application of A^^^ to a discrete sequence f[n] is concerned, 
the unfortunate difficulty is that the vectors bkU0 must necessarily lie on the 
lattice for (12) to be well-defined. In fact, it is easily seen that one cannot 
associate a digital filter with the RS operators in general. However, the good 
news is that, when N equals 2 and 4, the transformation f[n] ^ A^^ f[n] 
can be exactly realized without the need for interpolation by appropriately 
setting the scales bk of the directional running-sums. We will discuss the 
latter case in detail in §3. 

The algorithm for realizing (5) corresponding to a specified scale-vector 
map a(n) is based on relation (15). In particular, by considering the function 
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f{x) = f[n]S{x — n), and by proceeding exacdy along lines of the 1-D 
derivation, we express the filtered samples in (5) as 

2^-1 

f[n] = J2 ^i^i^ + r-Xi), (17) 

i=0 

where F{x) = ^ gh[n]f3^{x — n) denotes the interpolated version of the pre- 
integrated signal 5rb[n-] = Aj,^/[n]; r = 0.5{^{akin)-bk) cos Ok, J^iakin)- 
bk) sin^fe); and the pairs (xj, Wi) are the vertices and taps of the affine FD 
mesh in (16). Note that T,Wi and Xi are defined pointwise in (17); we 
dropped the index n to simplify the notation. We will discuss the implemen- 
tation aspects of the algorithm, particularly the computation of [n] and its 
interpolated form F{x) for the case N = Ain §3. 

2.2.3 Characterization of the radially-uniform box splines 

The motivation behind introducing the radially-uniform box splines was 
to develop eUiptical Gaussian-like filters, whose shape (size, elongation 
and orientation) can be continuously controlled, and a fast space-variant 
algorithm using such filters. Indeed, it turns out that the radially-uniform box 
splines (and its iterated versions) form close approximates of the Gaussian. 
To substantiate our claim we present the following result (proof in Appendix 
§6.1) that can well be seen as a "radial" version of the central limit theorem. 

Theorem 2.2. Consider the sequence of box splines /J^^g) i^)^ /^a(3) (^)' ■ ■ • 
responding to the scale-vectors a(2),a(3), . . where ak(N) = a^/{24/N) for 
1 < k < N. Then the following holds 




In fact, the radially-uniform box splines constructed using uniform scale- 
vectors are supported on a A^-sided uniform polygon, and it can be shown 
that they have continuous derivatives of order (A^ — 3). The above result 
is then consistent with the fact that the isotropy and smoothness of such 
box splines progressively improves with the increase in the directional-order. 
Moreover, it is also possible to mimic certain anisotropic Gaussians by using 
a sequence of non-uniform scale-vectors. Indeed, as a direct extension of 
Theorem 2.2, one can construct sequences of box splines which converge to 
anisotropic Gaussians as N increases. 



13 



Yet another useful form of anisotropic convergence is achievable based on 
the serial convolutions of a radially-uniform box spline, of a fixed directional- 
order, with itself. In particular, corresponding to fixed integers N and m 
(m > 1), and a scale-vector a = (ai, . . . , oat), we consider the iterated 
radially-uniform box spline 

fi^'^{x) = {^^*---*^^){x) (19) 

obtained through the (m — l)-fold convolution of P^{x) with itself. Then, 
for the particular sequence of box splines corresponding to the 

scale-vectors a(m) = (ai/-y/m, • • • , aN/\/rn), we have the result 



where 



1 

12 



(21) 



al cos^ Ok i flfe sin 29k 

Indeed, this follows directly from a certain version of the central limit theo- 
rem, which tells us that the each of the components 

converge to a "directional" Gaussian distribution as m — > oo. The covariance 
C is then given by the limiting sum of the covariances of the constituent 
box distributions. The utility of such iterated box splines will be discussed in 
§3.3. 

Having characterized the asymptotic behavior of the box splines, we now 
focus on the problem of approximating an anisotropic Gaussian using a box 
spline of a fixed directional-order. Since a centered Gaussian is uniquely spec- 
ified by its covariance, we propose a finite-order box spline approximation of 
the same based on its covariance. This, in fact, amounts to a moment-based 
approximation of the Gaussian, that is, the box spline resembles the Gaussian 
up to its second-order moments. Moreover, since the level-sets of Gaussians 
are ellipses, this equivalently amounts to constructing elliptical filters of 
different size, elongation, and orientation. 

The covariance of the radially-uniform box spline, namely. 



C^ = J xx'^P^{x)dx, 
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Figure 4: Intensity distribution of (a) the radially-uniform box spline, and 
(b) the separable B-spline, of order four. The respective scan profiles along 
7r/8 are shown in (c) and (d) . 
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(22) 



can be expressed as the sum of the covariances of the box distributions (cf. 
Appendix §6.2 for details) as follows: 

This provides the explicit dependence of on the scale-vector. In partic- 
ular, based on the eigen decomposition of C^, we propose the following 
characterization of the elliptical parameters of the box spline: Let Amax and 
Amin denote the largest and smallest eigenvalues of C^, and {vi,v2) the 
eigenvector corresponding to the eigenvalue Amax- The size s^, elongation 
, and orientation 9^ of the radially-uniform box spline fi^(x) are then 
defined as 



— Amax + Amin — _ „ ^ ^ Qfej 



= tan-^ ( - ) = tan-^ I ^ T\~Z' \ " ~ I , (23) 



12 

Amax _ S Cife + 
Amin Z]afc--/D' 



vij y J2al sin(26lfc 

where L>= {Y:alcos20kf + {Zalsm2ek)\ 

Since is strictly positive (see Appendix §6.2), all the above parameters 
are indeed well-defined. Note that the covariance matrix in (22) and the 
triple in (23) provide equivalent descriptions of the box spline geometry. The 
motivation behind introducing the latter is its convenient rotation-invariant 
nature: while changes with the rotations of a given box spline, and 

remains fixed. It is for this reason that we use the latter description for 
studying the dependence of the elliptical geometry on the scale-vector in the 
next section. 



3 FOUR-DIRECTIONAL BOX SPLINES 

We now study the particular /our-directionai box spline 

Piix) = (</?ai,0 * <Pa2,n/4 * </'a3,7r/2 * </^a4,37r/4)(^)' 

and the corresponding implementation aspects. This particular box spline is 
composed of patches of quadratic pol}aiomials (degree < 2), is continuously 
differentiable, and is compactly supported on a convex octagon (cf Fig. 3) . 
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{Xu, -w)^ ai Jxi5, +w) 




(Xl2, +w) 
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\{X7, —w) 

as 

){X3, +w) 



ixo,+w)^ OT 



Figure 5: The distribution of the taps of the FD mesh. The pairs {u, v) denote 
the positions {u) and the corresponding weights (v) of the taps of the FD 
mesh. 

We note that in [12] the authors have used separable B-spHnes to ap- 
proximate the Gaussian. Although these functions are built from the same 
constituent box distributions, the advantage of the four-directional box spline 
over the separable ones is that they are more isotropic. As seen in Fig. 4, the 
basic four-directional box spline, besides having a smaller support, exhibits a 
more Gaussian-like profile than the separable counterpart of identical order. 
In addition, the anisotropic four-directional box spline can be rotated to 
arbitrary orientations, while the separable ones are constrained to the image 
axes. 

3.1 Fast space-variant elliptical filtering 

The four-directional box spline is of particular interest in the context of the 
space-variant filtering following the fact that A^^^ can be implemented with- 
out interpolation when b = (1, \/2, 1, \/2). The corresponding interpolating 
function, (3^{x), turns out to be well-known in the box spline community, 
and is popularly referred to as the Zwart-Powell (ZP) element [2, 21]. The 
steps for realizing (5) using the four-directional box spline are as follows: 

(1) (Pre-integration) The crucial point is the choice of the scale-vector 
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b = (1, \/2, 1, \/2) corresponding to which the RS operator 

Ar*^ = ArJ, o A'i- , o Ar^ ,^ o A'i , 

b 1,0 V2,7r/4 1,t/2 V2,37r/4 

can be directly applied to f[n]. In particular, the running-sum gi,[n\ = 
Aj^^/[n] can be computed using the following four steps. 

(RSI) Horizontal running-sum: 

oo 

go[ni,n2] = A^^o/[ni,ra2] = ^f[ni - k,n2]. 

k=0 

(RS2) First-diagonal running-sum: 

oo 

g^/4[ni,n2] = A^^^^^c/o [m , "2] = V2^go[ni - k,n2 - k]. 

k=0 

(RS3) Vertical running-sum: 

00 

5^/2 [ni, 77,2] = A^^^/2f7r/4K,n2] = ^57r/4[m,n2 - k]. 

k=0 

(RS4) Second-diagonal running-sum: 

00 

9h[ni,n2] = ^^^37r/4f7r/2K,n2] = V2^g^/2[ni + k,n2-k]. 

k=0 

(2) (Finite-difference) At each position n, the FD mesh is computed using 
the scale-vector a(n). The weights Wi and the vertices Xi are listed in Table 1 
with the convention that a'^ = ak/V^ for k = 2,4. The mesh has a total of 
4x4 = 16 vertices; in particular, there are 4 clusters corresponding to the four 
boxes with 4 vertices per cluster, as shown in Fig. 5. The shift r = (n, r2) is 
given by Ti = (^/2ai+a2-a4-\/2)/2V2 andr2 = (a2+\/2a3+04-3\/2)/2A/2. 
The filtered sample is then computed using the formula 

15 

f[n]=J2^iFi^ + 'r-Xi). (24) 

i=0 

The interpolation samples F{x) = Y, 9b[n\l3^{x - n) in (24) are computed 
efficiently by taking advantage of the piecewise polynomial structure of the 
compactly supported ZP element (see Appendix 6.6). 
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Table 1 : Specification of the taps of the FD mesh associated with the operator 
A3. The weight w is given by (01020304)^^, where a = (ai, 02, 03, 04) is the 
corresponding scale-vector. 



i 




Wi 


i 




Wi 





(0,0) 


+w 


8 


(-04,4) 


—w 


1 


(ai,0) 


—w 


9 


(oi - O4, O4) 


+w 


2 


(a2>«2) 


—w 


10 


(o'2 -0^,0'2 + a^) 


+w 


3 


(01+02,02) 


+w 


11 


(oi + 03 — 04,02 + O4) 


—w 


4 


(0,03) 


— w 


12 


(-04,03 + 04) 


+w 


5 


(aijO's) 


+w 


13 


(01 - 04,03 + 04) 


—w 


6 


(02,03 + 02) 


+w 


14 


(O2 - 04, 03 + O2 + O4) 


—w 


7 


(01 + 02,03 + 02) 


—w 


15 


(01 + 02 — 04, 03 + o'a + 04) 


+w 



As in the 1-D setting, the running-sums are efficiently evaluated using 
recursions as summarized in Algorithm 1. The decisive computational advan- 
tage, especially for wider kernels, is derived from the fact that the number of 
vertices of the FD mesh is completely independent of the scale-vector. As a 
result, the algorithm has a fixed computational cost per pixel, modulo the 
cost of the running-sum and the interpolations (see Table 2). 

Algorithm 1 Space-variant elliptical filtering 

1. Input: /[n] and a(n) 

2. Perform recursions: 

9o[ni,n2] ^ f[ni,n2] +5o[m - ^,n2] 
9TT/4['ni,n2] ^ V^go[ni,n2] + 9n/A[ni - l,n2 - 1] 
9n/2[ni,n2] ^ 9^/4^1,1^2] + gn/2[ni,n2 - 1] 
9b[ni,n2] ^ V29T,/2[ni,n2] +5b[m + l,n2 - 1] 

3. Local FD operation: 
for each position n do 

compute Wi, Xi and r using a(n) 
evaluate F{n + t — Xi) using ZP interpolation 
/W ^ Ei WiF{n + T-Xi) 
end for 

4. Return f[n\ 
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3.2 Size, elongation and orientation of the box splines 

As was mentioned earlier, the size and shape of the radially-uniform box 
spline can be controlled by appropriately adjusting the scales of the con- 
stituent box distributions. In this regard, we now discuss the following: (i) 
the forward problem of controlling the anisotropy of the four directional 
box spline by acting on the scale-vector, and (ii) the inverse problem of 
uniquely specifying the scale-vector of the box spline corresponding to a 
given covariance (geometry). For notational ease, we shall henceforth drop 
the superscript = 4 when referring to the four-directional box spline and 
its related parameters. 



3.2.1 Control on the anisotropy 

The elliptical geometry of this box spline is specified using parameters defined 
in (23), namely, 

^a-^2^a„ 0a-tan I I , and - 



where D = (a| — a\Y + (^2 ~ ^I)^- turns out that the size and orientation 
can be arbitrarily controlled by adjusting the scale-vector. Indeed, the size 
can be easily manipulated by multiplying the scale-vector a with a constant 
factor, since this leaves both the orientation and elongation unchanged. The 
elongation can be arbitrarily controlled in the neighborhood of the four 
principal directions. However, there exists a finite upper bound on the 
elongation along other directions (cf. Appendix §6.3). 

Proposition 3.1. For every (f> in [0, tt), there exists a scale-vector a such that 
9 Si = (j). There is however a finite bound on the elongation, and is given by 

1 + + 

sup ^a < U{4>) = 1 , (25) 

where = ^(tan (p — cot (^)sign(| — cf). The supremum is over the set of a for 
which 9a = (j)- 

Fig. 6 illustrates the variation of l/U{4)) as a function of the orienta- 
tion; the rationale behind showing the inverse plot is to avoid the blowups 

U{4>) — > +00 as (j) — > 9k. In particular, a bound of 3 + 2\/2 ^ 5.8 is attained 
along the orientations (p = {2k — l)7r/8, 1 < fc < 4, exactly mid-way between 
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Figure 6: Polar plot of the symmetric variation of !/[/((/)) as a function of the 
filter orientation 0, where J7(</>) is the bound on the elongation. The bound 
reaches its minimum when the orientation of the filter is exactly midway 
between two principal axes, whereas arbitrary elongation is achievable in 
the neighborhood of the four principal directions (/> = 0, 7r/4, 7r/2 and 3tt/A. 

two adjacent primal directions. This is perfectly reasonable since the control 
on the geometry of the box spline is minimal along these directions. 

In order to specify the elliptical geometry of the box spline, we use either 
of the following equivalent descriptors as per convenience: 

(Dl) Size, elongation and orientation {s, g,6). 

(D2) Length of the major and minor axes, and the orientation (A/Amax, V^min, d)- 
(D3) Covariance matrix C. 

Descriptor (Dl) stipulates the lengths of the major and minor axes as 
Amax = sg/{l + q) and Amm = + g), respectively, whereas, (D2) gives 
the corresponding covariance as 

Q ^ ( -^max COS^ 9 + Amin sin^ 6 \ (Amax " Amin) sin 26 \ 
\ 5 (Amax - Amin) sin 26* Amin COS^ + Amax sin^ ]' 
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3.2.2 Optimal scale-vector for a given anisotropy 

Since the covariance matrix of the four-directional box spline (c£ Eqn. (22)) 
is given by 

^ 24 V al-a\ 2al + al + al J ' ^ ^ 

the inverse problem is that of specifying a scale-vector a such that Ca = C. 
By introducing the positive vector p = (af , a^, a|, 04), the problem can be 
reformulated as: find p > 0, such that Mp = c, where 

/ 2 1 1 \ 

M= 10-1, and c = 24(c(l, 1), C(l, 2), C(2, 2)) . 
\ 1 2 1 / 

The scale-vector solution is then given by a = y^. As far as existence 
of solutions is concerned, proposition 3.1 ensures that the linear system 

Mp = c, p > 0, corresponding to a given geometry (Amin, \nax, 0), is always 
solvable provided that the elongation q < U{6). Moreover, as it turns out, 
the linear system is under-determined and has infinitely many solutions. The 
idea then would be to use a scale-vector that is "optimal" in some sense. But 
first, we try to characterize the solution space of the system Mp = c, p > el. 
For reasons that will soon be obvious, we propose to modify the positivity 
constraint as p > el, where e is some arbitrarily small positive number. We 
observe that M is full-rank, and hence the null-space is of dimension 4— 3 = 1. 
In particular, this signifies that the solutions of Mp = c lie on the affine 
subspace {p + te : t G R}, where p is a particular solution (Mp = c) and e is 
in the null-space (Me = 0) . Moreover, one can easily verify that in order to 
adhere to the positivity constraint, p+te = {Pi+t,p2—t,p3+t,p4—t) > el, it 
is both necessary and sufficient that t lies in the closed interval [ti, tr], where 
= max(— pi + e, —ps + e) and U = min(p2 — e,p4 — e). In particular, we set 
e = (1, —1,1,-1) which is in the null-space of M. Note that one can easily 
compute p by pivoting one of its components and solving for the remaining 
three; since M is of full-rank, the reduced system is always solvable. We now 
use the available degree of freedom to select a solution that maximizes a 
certain measure of Gaussianity. 

A classical measure of the Gaussianity of a 1-D function is its kurtosis 
(the fourth-order cumulant). For a centered function /(x), this is defined 
as K = ^4 — 3fi2, where fi4 and fi2 are the fourth-order and second-order 
moments of f{x), respectively. The central property of this measure is that 
K = for a true Gaussian function, and as a result, the absolute value of 
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(a) 9 = (b) 6 = 7r/9 (c) 6 = 37r/4 (d) 6 = 7n/18 




{e)0 = n/2 (/) f = 237r/36 {g)9 = 3Tv/A (/i) 6* = 87r/9 



Figure 7: Intensity distributions of the four-directional box splines of identical 
size (s = 1) and elongation (g = 2.5), but with different orientations. The 
ellipse in each figure represents a level-set of the Gaussian having the same 
covariance as the corresponding box spline. 



23 



the kurtosis provides a measure of Gaussianity of the function. In particular, 
smaller absolute values correspond to more Gaussian-like functions. 

As for a bivariate function f{x), we shall use the following matrix-valued 
extension 

K = L - tr(C)C - 2C^ (27) 

where C = j{xx^)f{x)dx and L = J {xx'^Y f{x)dx are the second-order 
and fourth-order moment matrices of f{x), respectively [17]. This consti- 
tutes a valid extension of the 1-D kurtosis since (27) reduces to k = /7,4 — S/x^ 
when d=l. Moreover, we also have the following desirable properties: 

(i) If f{x) is a multivariate Gaussian, then K = (cf. [17] for a proof). 

(ii) The Frobenius norm of K, namely, 

l|K|| = (E|K(z,,)|^)''' 

is rotation-invariant, i.e., the kurtosis matrices of the rotations of /(x) have 
the same Frobenius norm (proof in §6.4). 

Following the above arguments, we propose to solve the optimization 
problem 

Po = argminp ||Kp||^, Mp = c, p > el. (28) 

This jaelds the optimal scale-vector ao = ^/po corresponding to the most 
Gaussian-like box spline. The rotation-invariance property ensures that the 
box splines of identical size and elongation but different orientation, obtained 
via the solutions of the above optimization problem, are as homogenous as 
possible. 

The norm of the kurtosis matrix of I3a{x) turns out to be ||Kp|p = 
'^kPk + + pI){P2 + pI) (see Appendix §6.5). Substituting Pk=Pk + teu 
into this expression, we arrive at the quartic polynomial 

CW = Y^^Pk + ektf + {{pi+ tf + (P3 + tf) {{P2 - tf + (P4 - tf] , 
k 

which, together with the parameterization p = p + te, simplifies the problem 
to one of finding 

to = arg minj ({t), t e [te, U]. (29) 

The optimal solution is then given by ao = y/p + toe. This problem however 
is easily solved, since the minimum is attained either at one of the interior 
points (te, U) where ('{t) = 0, or at one of the boundary points. In particular. 
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we have the following simple algorithm for designing optimized Gaussian- 
like box splines of a specified covariance: 

(i) Set p4 = 1, and compute p by solving the sytem Mp = c. 

(ii) Use p to compute t^, U and the coefficients of C'{t). 

(iii) Find the real roots of = over the interval {ti,tr); denote the set 
of real roots by R. Then ao = (p + toe)^^'^, where^ to = argmin^ C(*))^ € 
RU{te,tr}. 

In particular, the coefficients of the cubic equation ('{t) = Cit^ + (2^^ + 
Cst + C4 = in (iii) are given by Ci = 32, C2 = 24(pi - P2 + Ps - P4), Cs = 
16 EPfc - 8(pi +P3){P2 +P4), C4 = 4(pf -pI +pI -pI) + 2(pi +p3){pl +pI) - 
'2{p2+p4){pi+Pl). 

The box splines obtained using the above optimization at various orienta- 
tions are shown in Fig. 7. The quality of the Gaussian approximation under 
different practical settings of the orientation and the elongation is quantified 
in Fig. 8. 

The correspondences (1, g, 9) o (ai, 02, as, 04) for < < vr and I < g < 
U{6) can be pre-computed and stored in a look-up table. Note that for a given 
g, the set of correspondences (1, g, 6) o (ai, 02, as, 04) have an inherent four- 
fold symmetry in 6 owing to the presence of the four principal directions. 
Hence, it suffices to store the scale-vector correspondences for < < 7r/4 
which reduces the size of the LUT by a factor of four Indeed, for any 
arbitrary size s > 1, orientations < ^ < tt, and elongation 1 < ^ < U{d), the 
corresponding scale-vector is then obtained through the following operations: 

(01) Rotation: 

'e for 0<e<7r/4 

,_\o-tt/A for 7r/4 < < tt/2 
I ^ - 7r/2 for 7r/2 < ^ < 37r/4 
^^-37r/4 for 37r/4 < ^ < tt. 

(02) Find (ai, 02, 03, 04) corresponding to (1, g, cp) using the LUT. The desired 

^The tie is randomly broken if ({t) has multiple minimizers over [ti,tr\ (this was rarely 
reported in practice). 
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Figure 8: Normalized correlation between the optimal four-directional box 
spline and the target Gaussian at different elongations and orientations. For 
a fixed elongation, the correlation is minimum at the critical orientation 6 = 
22.5°, and improves symmetrically as 9 approaches the principal orientations 
(cf. Fig. 6). 

scale-vector is then given by the following permutation and rescaling: 



(01,02,03,04) I-)- < 



^(01,02,03,04) 
^(02,03,04,01) 
^(03,04,01,02) 
^(04,01,02,03) 



for O<0 < 7r/4 

for TT/A<e < Tr/2 

for TT/2<e < 37r/4 

for 37r/4 <0 <tt. 



3.3 Higher-order box splines 

As suggested by the convergence result (18), the Gaussian-like nature of 
the four-directional box splines can be improved by using more directions. 
Implementing the corresponding space-variant filtering using the algorithm 
in §2.2.2 however turns out to be challenging and not very practical — the 
principal axes of these box splines are generally along off-grid directions, 
and one needs to interpolate the image for implementing the associated 
running-sums. 

The iterated four-directional box splines /3a™(a3) introduced in §2.2.3 
provide a practical alternative. These box splines rapidly converge to a 
Gaussian with the increase in m. Also, note that the four-directional box 
spline and its iterates have identical covariances. This implies that the 
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(a) Isotropic forms (s = 1, ^ = 1) 



(b) Anisotropic forms (s = 1, g = 3,9 = tt/S) 

Figure 9: Higher-order box splines through iterative convolutions. Left: The 
reference four-directional box spline; Center: Iterated box spline obtained 
by convolving the (reseated) four-directional box spline with itself; Right: 
Target Gaussian having identical covariance. 
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algorithm in §3.2.2 can be used for optimizing the iterated box sphnes as 
well. The first two iterates of the four-directional box spline along with the 
target Gaussian are shown in Fig. 9. It is seen that (3a'^{x) resembles the 
target Gaussian very closely. In fact, the minimum correlation coefficient rises 
from 95% to 99% for m = 2 (cf Fig. 8). In practice, we can thus implement 
a higher-order Gaussian-like filtering by simple iterations of the algorithm in 
§3.1. It suffices to set the scale-vector in the algorithm as a./^/rn, where m is 
the number of iterations. 

4 EXPERIMENTAL RESULTS 
4.1 Computation time 

The space-variant filtering using the four-directional box spline was imple- 
mented in Java on a 2.66 GHz Intel system. The typical execution times 
required for convolving a 512 x 512 image with kernels of various sizes are 
shown in Table 2. It is clear that the run time is independent of the size of 
the kernel. 

Table 2: Average computation time for box splines of different sizes. 



Size (s) 


1 2 4 8 16 


Time (millisec.) 


101 100 103 101 100 



4.2 Application: Feature-preserving smoothing 

We now present an application to demonstrate the space-variant algorithm 

described in §3.1. Filtering of noisy images using isotropic Gaussian filters 
often results in excessive blurring of the anisotropic image features. Diffusion 
filters are known to perform better in such cases [19]. As an alternative, 
we propose to filter the corrupted image using our anisotropic Gaussian-like 
filters, where we adapt the size, elongation and orientation of the filter to 
the local image features. The main idea is to locally average the image 
using elliptical windows that have been elongated along the image feature 
(orthogonal to the local gradient). This induces more smoothing along the 
direction of minimal intensity variation resulting in the suppression of the 
ambient noise, while preserving the sharpness of the image features. 

To derive an estimate of the local image anisotropy, we use the paradigm 
oi structure tensors [9], where the local orientation 6{x) is estimated through 
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the minimization of a certain weighted norm of the directional derivative. 
In particular, if we denote the directional derivative of f{x) along along 
U0 = (cos 6, sin 0) by Dof{x), then 8{x) is given by the minimizer of 



/ w{s)\{Def){x-s)\^ds 



(30) 



where Q is the support of the isotropic averaging window w{x). The solubility 
of the above optimization problem follows from the observation that this 
can be recast as an eigenvalue problem. In particular, by expressing the 
directional derivative in terms of the gradient g(£c), namely as DQf{x) = 
UQg{x), one can rewrite (30) as 



where the structure tensor J (a;) is the 2x2 positive-definite matrix 



The local orientation 0*{x) is then given by the minimizer of (31) associated 
with the minimum eigenvalue of J (a;). 

In view of the definitions in (23) and the fact that the eigenvalues of J (a;) 
are always non-negative, we propose to estimate the elongation as follows: 
We set Q* = Amax/Amin if both eigenvalues are non-zero, equal to 1 if both 
are zero (locally isotropic intensity), and equal to max(l, A) if only one of the 
eigenvalues A is non-zero. Finally, we estimate the size of the box spline as 
s* = Amax + Amin- The triple (s*, Q*,9*) is then used to compute the optimal 
scale-vector using the algorithm described in §3.2. The components of J can 
be efficiently computed using simple convolution and pointwise operations; 
we refer the reader to [9, Chapter 13] for implementation details. The main 
steps of the proposed smoothing algorithm are: 

• Computation of the structure-tensor. 

• Pre-integration of the corrupted image using the running-sum filters. 

• Computation of the triple {s*, g*,9*) at every feature location using 
the structure tensor This is used to compute the scale-vector a(n) 
of the optimal Gaussian-like box spline using the algorithm in §3.2. 
Isotropic box splines are used in the uniform-intensity regions; we set 
a(n) = {a, a, a, a), where a is proportional to the noise variance. 



UQ3{x)ue 



(31) 
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Figure 10: Results on a test image. (A) Barbara corrupted with additive Gaus- 
sian noise, PSNR = 18.0 dB; (B) Isotropic smoothing, PSNR = 23.10 dB; (C) 
Diffusion filtering, PSNR = 23.25 dB; (D) Our algorithm, PSNR = 23.58 dB. 
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(C) (D) (E) 




Zoom (C) Zoom (D) Zoom (E) 



Figure 11: Results on a real image. (A) Noise-free immunofluorescence 
image of actin fibres (Courtesy of C. Aemisegger, CMIA, University of Zrich); 
(B) Image corrupted with additive Gaussian noise, PSNR = 12.20 dB; (C) 
Isotropic smoothing, PSNR = 15.38 dB; (D) Diffusion filtering, PSNR = 
15.50 dB; (E) Our algorithm, PSNR = 15.80 dB. 
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• Computation of the FD mesh using a(n), and its appUcation to the 
pre-integrated image. 

To demonstrate the effectiveness of our strategy in preserving oriented 
patterns in noisy images, we compare the results obtained from our algorithm 
with those obtained using the (fixed-scale) isotropic Gaussian filter and the 
Perona-Malik diffusion filter [14]. We use the standard test image oi Barbara 
and corrupt it with additive Gaussian noise. The variance of the noise is used 
to set the size of the Gaussian for the isotropic smoothing. The parameters 
used for the Perona-Malik filter were typical: time step of 0.1, conductance in 
the range of 10 30, and a total of 15 ~ 30 iterations. The parameters were 
manually tuned to optimize the PSNR, and also to avoid blocking artifacts. 
Fig. 10 shows the results obtained from the different filters. As far as the 
quantitative evaluation of the filters is concerned, our algorithm clearly 
outperforms both isotropic and diffusion filters in terms of the Peak-Signal- 
to-Noise-Ratio (PSNR). Moreover, as shown in the zoomed-in sections of 
the respective images, the oriented stripes on the clothes are quite faithfully 
restored by our algorithm. A significant amount of blurring of the stripes 
is seen in the results obtained using isotropic and diffusion filtering. The 
non-linear diffusion filter, however, tends to perform better at low PSNRs in 
the range of 5-10 dB (cf Table 3). 

Next, we compare the results on a real biological image and at a much 
lower PSNR of around 12 dB. We consider the fluorescence image shown 
in Fig. 11, which exhibits numerous elongated fiber-like structures. The 
parameters of the isotropic filter and the diffusion filter are set as in the 
previous case, except that the iteration count for the latter is increased to 15. 
As before, the improvement of the PSNR obtained using our filter is higher. 
Importantly, as seen from the zooms, our algorithm results in significantly 
less merging of the close fibers and blurring of the finer ones. The average 
execution time of our algorithm is 0.6 seconds for a 512 x 512 image, which 
includes the computation of the structure-tensor, the running-sums, the 
optimal scale-vector, the interpolated samples and the finite- differences. 

The four-directional box splines can also be used to derive fast space- 
variant detectors based on Gaussian forms, e.g., the Laplacian-of-the-Gaussian 
(LoG) or the so-called Mexican-hat detector. We refer the interested reader to 
[3], where the isotropic forms of the four-directional box spline were used to 
realize a fast and scalable Mexican-hat-like detector. In particular, a modified 
version of the space-variant algorithm described in §3.1 is used to design 
an efficient coarse-to-fine strategy for the detection of centers and radii of 
cells/nuclei in fluorescence images. 
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Table 3: Comparison of the filters at different noise levels using the test 
image of Barbara. The table shows the PSNR of the outputs. 



Input PSNR (dB) 


10.0 12.0 14.0 16.0 18.0 20.0 


Isotropic filter 


15.38 18.20 20.20 21.65 23.10 24.30 


Diffusion filter 


15.48 18.31 20.30 21.70 23.25 24.35 


Our filter 


15.45 18.38 20.57 21.94 23.58 24.56 



5 CONCLUSION 

In this paper, we presented a framework for elliptical filtering using the 
radially-uniform box splines. The associated space-variant filtering was 
efficiently realized using running-sums and local finite-differences. The 
attractive features of our algorithm are: 

• the 0(1) computational complexity per pixel, and 

• the use of real-valued parameters for continuously controlling the shape 
and size of the filter. 

Our filtering paradigm offers a nice trade-off between the quality of approxi- 
mation of Gaussians and the computational complexity of linear space-variant 
filtering. 

We also presented a closed form solution for the problem of constructing 
four-directional box splines with given covariances. The scope of our algo- 
rithm was demonstrated through the realization of a smoothing filter that 
can adapt to the local image characteristics. 

6 APPENDIX 

6.1 Proof of theorem 2.2 

We first establish that the Fourier sequence /3a(2) "^(3) ('*')'■■ • conver- 
gences pointwise to a Gaussian: 

^lim^^^^)M = exp(-^^||c.f). (32) 

We then show that the above convergence is also in the L^(R^) norm. This 
will establish the theorem, since it is well-known that the Fourier transform 
of a Gaussian is a Gaussian, and that — y g in if Jn — > y in L^. 
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To derive (32), we note that ^a,e('*^) = PaiugU}) = sine [augu/2), 
where siiic(x) = sin(x)/x for x ^ 0, and equals 1 at the origin. Then, 
the convolution-multipHcation rule gives 

i^) = n <me, i^) = U sine ("^ulu) . (33) 
k=i fc=i V ^ / 

Using the estimate sinc(x) = 1 — + O(x^) for \x\ < 1, and substituting 
afc(iV) = a^/{2A/N) into (33), we have 

/ 2 \ 

^2iV)(^) = n 1 - ^K'^)' + O {N-') {M < cN) (34) 

k=l ^ ^ 

where c is some positive constant. By developing the quadratic factor (u^^ijS)^ 
and the product in (34), we arrive at the estimate 

^( 2 2 2 ^ 

! 1 ^ 



k=\ 

1_^||^||2) +^(^2_ 2)^_^||^||2\ y-cOs2^, 

2iV " " / 2iV^ ^ V 2iV " " / ^ " 

2 / 2 \ ^-1 ^ 

+ ^'"^'^^ ^ " ^ ii'^'in ^ + ^ 
^ fe=i 

= (1 - 1^ ll'^f + O (iV-2) (||u;|| < cN). (35) 

This is exactly where the fact that are uniformly distributed over [0, vr) is 
invoked: the cancellation of the linear factors in the second step is based 
on the identities Y^=\ cos 20^ = 0, and Y^=\ sin 2^fc = 0, where = 
(fc ^ l)7r/A^. Since (1 ~x/m)"^ converges to cxp(— .t) as m — )■ 00, it can now 
be readily seen that (32) follows as the limiting case of (35) . 

To demonstrate that (32) holds in the norm sense, it suffices to 
show the sequence of error functions £^Ar(a;) = (i^^^^^.^iu) — exp(— cr^ ||u;||^ /2) 
converge to zero in the above norm, i.e., ||£^Ar||L2 — > Oas N — > 00. Since we 
have already demonstrated that £n{<^) — > pointwise, all we need to show 
in order to invoke the dominated convergence theorem is that the sequence 
|f2('*')| , |^3('*')| , • • • is uniformly bounded by a function. Moreover, since 

|^iv(u;)| < +exp(-a2 ||a;f /2), 
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it, in fact, suffices to show that each admits such a bound. 

The main idea behind establishing such a bound is that the above men- 
tioned sequence can be covered by a Gaussian in a neighborhood of the 
origin and by a function with sufficient decay at the tails, both of which 
are independent of N. Indeed, using the estimate sinc('u) < 1 — u^/tt^ for 
|u| < TT, one can verify that 



N 



k=l 

N 



6a2 



< 



n 1 



k=l 



T 

ffc 



) 



<exp(-Ci||a;f ) {\M < S) 
As far as the tail is concerned, the Cauchy-Schwarz inequality, < 



l^^fcll • = gives 

N 



k=l 



< 



Co 



a; 



|a;|| > S). 



Here Ci, C2 and 6 are appropriate positive constants that are independent of 
N. Combining the above estimates, we see that 



|/3^^)M|<exp(-Ci||u»r) + 



Co 



\u:\\ 



1 — rect 



\U}\ 



for all a;. Since the function on the right is indeed in L^(R^), this establishes 
the desired bound, and consequently, the norm convergence. 



6.2 Covariance matrix 

We begin with the observation that if f{x) and g{x) are symmetric (about 
the origin) and have a total mass of unity, then C/*^ = C/ + Cg, where C / 
denotes the covariance matrix of f{x). Indeed, by noting that /(O) = g{0) = 
1 (unit mass) and 5^/(0) = ^^^(O) = (by symmetry), and by recaUing the 
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multiplication-differentiation rule / XiXjf{x)dx = —didjf{0), we have 
Cf*gii,j) = J XiXj{f * g){x)dx 

= -mdidjfio) - mdid^m - dimdjm - dimdjm 

= -didjfio) - didjm 

= Cf{i,j) + Cg{i,j). 

Since the directional box distributions ipa^fi^{x) satisfy these criteria, we 
have that = C^, where is the covariance matrix of ipa^fi^i'^)- We 
explicitly compute the component C(l, 2); the remaining components can 
be similarly derived. Using the multiplication-differentiation rule again, we 
have 



Cfc(l,2) = j xiX2(pak,eki^)dx 



7? 

sin 29k. 



w=0 



24 

Therefore, (1,2) = EfeCfc(l,2) = ^fc 4 sm 2^^/24. 

As far as the positive-definite nature of is concerned, it will suffice to 
show that its eigenvalues 

Amax =\{^4 + ^) and Aniin = ^ 4 " ^) 

where D = {J^a-l cos(26'fc))^ + {J2 al sin(20fc))2, are strictly positive. This is 
obviously the case for Amax- Moreover, the inequality 

( E 4) ' - = ( E ' - ( E ' - ( E 4 sin(20.)) ' 

k k k k 

= 2 ^ alaj (l - cos(26'fe - 26*^)) > 

k^£ 

tells us that X) > Vd. Hence, Amin is strictly positive as well. 



6.3 Proof of proposition 3.1 

Following definition (23), the dependence of orientation of the box spline 
/3^{x) on the scale-vector can be expressed as 

tan ea = 1^ + sign(a2 - 04) V^l+~z^ {0 < 6 < tt) (36) 
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where i/ = (a| — al)/{al — al). 

We note the following: it is both necessary and sufficient that 02 > 04 
(resp. 02 < 04) for the box spline to be oriented between < 6*3 < vr/2 
(resp. 7r/2 < < tt), and it is the map (01,02,03,04) (z^, sign(a2 — 04)) 
that uniquely determines the orientation of the box spline. Indeed, the 
uniqueness aspect is based on the argument that, for < < tt/2, (36) 
reduces to tan 9^ = 1^ + vTT^ as a consequence of the necessary condition 
02 > 04. This implicitly represents a one-to-one between 0^ and u over the 
domains (0, 7r/2) and (—00, 00), since the map 9^ i-> tan6'a from (0, tt/2) into 
(0, 00), and the map u 1-^ 1/ + \/l + from (—00, 00) into (0, 00) are both 
strictly monotonic. In a similar vein, a one-to-one between and u over the 
domains {tt/2, tt) and (—00, 00) can be established. In particular, we have 

v = ^(tan6'a - cot 6'a)sign ^ ^a) • (37) 

That is, given any orientation 9^ = (p, the corresponding is uniquely 
specified by (37). This establishes the first part of the proposition, since there 
trivially exists some positive vector (oi, . . . , 04) such that (a|— of )/(a|— 04) = 

As far as the bound is concerned, we observe that the elongation can be 
expressed as ^a = 1 + 2/(7 — 1), where 7 = X)o|/\/l)a > 1- For a given 
orientation 9^ = </>, the components of the feasible scale-vectors bear the 
relation (of — a|) = i'^{al — a|), and thus we have that 



7 



E4 



al + al ^ 03 + 04 



1 of + oj ^ \u<p\ al + al ^ 1 + W4,\ 



following the trivial inequalities of + o| > |of — a||, and al + al> {a^ — al\. 
Consequently, 

= 1 + T < ^, (38) 

7-1 1 + 1^^1- /TT^ 



The above bound is tight since it can be approached arbitrary closely by 
making the scales o^ and ak {9i < 4> < 9k) arbitrarily large. 
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6.4 Rotation-invariance 

Let K and K^i denote the kurtosis matrices of /(cc) and its rotation /(Hqx), 
respectively, where is the rotation matrix. Observe that the matrices Lo 
and L are related as 



Similarly, we have Ce = R^CRj. This also gives us the equivalence tv{Cg) = 
tr(ReCR^) = ti^CRjKe) = tr(C) following the identities tr(AB) = 
tr(BA) and RjRe = I. We can then write 

Ke = L0- ti{Ce)Ce - 2C^ReLRj - tr(C)R0CRj - 2ReC2R^ 
= ReCL - tr(C)C - 2C2)R^ = R^KR^. 

Our claim follows immediately, since ||K5(||^ = tr(K^K5i) = tr(RgK^KRg') = 
tr(K^K) = ||Kf . 

6.5 Kurtosis measure 

In order to compute the kurtosis matrix, we only need to evaluate the fourth- 
order moments; the second-order moments are already known. In particular, 
using Fourier identities similar to the ones used in §6.2, one can derive the 
following expressions: 

xf/3s,{x)dx = -fj,4^ [iaf + 02 + a|) + [Gafal + dalal + 3a|a|) , 

/I 3 
xlx2Pa{x)dx = -/X4 (02 " "I) + 21^2^1 (^2 - af) > 

/Q. Q, rt / \ 1 /A 4\ "^2/22 22 22 S2 22 ^2 2\ 

XiX2Pa{X)dX = -fX4 + a^j ['^1^2 + ^l'^4 + (^2^3 + ^3'^4 ~ '^2'^4 + ^(^l^^a) > 

/I 3 
xixl/3aix)dx = -H4 (^2 - at) + 21^1(4 (^2 - al) > 

j X2/3a{x)dx = ^/X4 {iaj + 02 + af) + {Qalaj + 6a|o4 + 30304) , 

where ^4 = 1/80 and ^2 = 1/12 denote the fourth and second-order mo- 
ments of respectively. These provide the components of the matrix La, 




{y = Rix) 



= R^iLR^. 
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which in turn leads to the following simple expression for the kurtosis matrix 



Ka = La-tr(Ca)Ca-2C2 



( 



(39) 



Finally, from (39), we get ||Kaf = ELi 4 + {af + 4)(4 + aj). 

We note that the negative factor (/i4 — Sfi^) in (39) is in fact the kurtosis 
of rect(x), the sub-Gaussian constituent of the box spline. The fact that Ka 
is negative-definite is thus consistent with the sub-Gaussian nature of the 
resulting box spline. 

6.6 Fast ZP interpolation 

Given a discrete function c[n] and a point x on R^, we outline a technique 
for the fast evaluation of the sum 



where b = (1, ^/2, 1, ^/2). A sketch of the partitions of the piecewise polyno- 
mial (3^{x) is provided in Fig. 12. The exact functional forms of the ZP box 
spline 2/3^(a;) corresponding to these partitions can be found in [21]. Since 
/3^(a;) has a compact support, this is in fact a finite sum, and requires at most 
seven evaluations of the function /3^(- — x) for any arbitrary translation x. 
This is illustrated in Fig. 12, where the red dots xo,. . . ,xq denote the lattice 
points that intersect the support of f3^{- — x). Thus, one needs to evaluate 
the translated ZP at the points xq, ■ ■ ■ ,xg in order to compute the sum. The 
drawback here is that naive evaluation of the spline at every xj requires a 
decision-making to figure out the associated partition before computing the 
corresponding polynomial. 

The redundancy that we exploit is as follows: Consider the triangular 
regions Pq, . . . ,P3 marked with blue dashed lines in Fig. 12 corresponding to 
the four different partitions of the ZP. These together constitute a unit cell of 
the lattice, and hence only one lattice point intersects them. The figure shows 
a particular instance where this point, denoted by xq, lies in Pq. This clearly 
fixes the partitions of the remaining lattice points xi, . . . , xq. Thus, if we de- 
note the polynomials corresponding to these partitions by pQfi{x), . . . , pq^q{x), 
then the sum in (40) is simply given by Yl%o '^[^j]Po,ji^j)- More generally, if 
Xq intersects the partition Pi {0 < i < 3), and if we denote the corresponding 
polynomials by pij{x), then the sum is given by Yl^=o'^[^j]Pi,j{^j)- Thus, 
we have the computational advantage that at most two binary decisions 




(40) 
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Figure 12: The figure shows the translated ZP box spHne j3^{- + x),h = 
(1, a/2, 1, a/2). The red dots xi, . . . ,xj correspond to the points on the Carte- 
sian lattice, and the triangular regions Pi , . . . , P4 are different partitions of 
the ZP, which together constitute a unit cell of the lattice. 

are required to simultaneously determine the ZP partitions corresponding 
to the points Xj, where the coefficients of the polynomials Pij{x) can be 
pre-computed. 
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